Non-axisymmetric oscillations of rapidly rotating relativistic stars by conformal 

flatness approximation 
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We present a new numerical code to compute non-axisymmetric eigenmodes ol rapidly rotating 
relativistic stars by adopting spatially conformally flat approximation of general relativity. The ap- 
proximation suppresses the radiative degree of freedom of relativistic gravity and the field equations 
are cast into a set of elliptic equations. The code is tested against the low-order f- and p-modes of 
slowly rotating stars for which a good agreement is observed in frequencies computed by our new 
code and those computed by the full theory. Entire sequences of the low order counter-rotating 
f-modes are computed, which are susceptible to an instability driven by gravitational radiation. 



INTRODUCTION 

Computing characteristic oscillations of rapidly rotat- 
ing stars in general relativity (GR) is much harder than 
in Newtonian theory where it is already a difficult task 
because of deformations of equilibria from a spherical fig- 
ure. Relativistic gravity is expressed by a spacetime met- 
ric tensor which has its own dynamical degrees of free- 
dom, compared to a single scalar potential in Newtonian 
gravity. These dynamical degrees of freedom are encoded 
in the theory in such an intricate way that they may be 
separated from non-dynamical ones only in cases with 
a high degree of symmetry. The equilibrium far from a 
sphere raises another important issue of how to impose 
boundary conditions for gravitational radiation. Free os- 
cillations of a star in GR are characterized by a condition 
of gravitational waves that propagate "purely outward" 
at the infinity. For rapidly rotating stars, there have been 
so far no prescription to impose this condition. 

Up to now the only linear eigenmode problem of 
rapidly rotating stars solved in full GR is that of the 
neutral mode when the eigenfrequency vanishes (Ster- 
gioulas and Friedman [H, hereafter SF98). For rotat- 
ing stars in post-Newtonian theory Cutler and Lind- 
blom performed eigenmode analysis and computed 
sequences of low order f-modes. Apart from them the 
traditional eigenmode problem of oscillations of rapidly 
rotating stars had been solved only within the Cowling 
approximation 0] , where Eulerian perturbations of met- 
ric coefficients are neglected. 

Remarkable progresses in numerical relativistic hydro- 
dynamics, however, have made direct simulations of stel- 
lar oscillations possible. The Thessaloniki-MPA group 
computed oscillations in the Cowling approximation! 4] 
and in the spatially conformally flat approximation [5j. 
The AEI group computed non-linear oscillations of rotat- 
ing stars in full GR [6]. The Tubingen group has devel- 
oped linear evolution codes specialized for oscillations of 
rotating stars with the Cowling approximation Q. The 
Tiibingcn-LSU-Thcssaloniki group computed the low or- 



der f-mode sequences in full GR [8j. 

Despite the remarkable successes in computing oscilla- 
tions of relativistic rotating stars by direct hydrodynam- 
ics simulations, the traditional eigenmode problem still 
needs to be studied. Firstly extracting eigenmodes by 
dynamical simulations needs much more computational 
resources than the eigenmode analysis. Especially the 
extraction of (the modulus of) an eigenfunction needs 
multiple numerical runs to obtain a single mode [4] . Sec- 
ondly it is not straightforward to identify eigenmodes by 
looking at the results of dynamical simulations. As the 
observational asteroseismology [§| the knowledge of the 
eigenmode analysis is indispensable to interpret what is 
obtained by a numerical experiment. 

In this paper we present a non-axisymmetric eigen- 
mode analysis of rotating relativistic stars beyond the 
Cowling approximation. For eigenmodes with low order 
and degree, which may be mainly excited in an astrophys- 
ical situation and are of interest in gravitational wave as- 
tronomy, the Cowling approximation works poorly with 
large errors Thus we here take into account gravita- 
tional perturbation with suppressing gravitational radia- 
tion. For typical compactness of a neutron star (~ 0.2), 
the damping time of a typical fluid oscillation due to 
gravitational radiation is three orders of magnitude larger 
than the oscillation period (e.g., Andersson and Kokkotas 
[l(|) and it may be neglected in evaluating the charac- 
teristic oscillation frequencies. 



FORMULATION 

In the context of relativistic astrophysics, approxima- 
tions of GR with the omission gravitational waves are 
considered by Isenberg (Isenberg [llj ]. originally writ- 
ten in 1978). The idea is that the behavior of a grav- 
itating system may be well-approximated without grav- 
itational emission as far as its timescale is longer than 
the dynamical one. Isenberg proposed several recipes 
to eliminate the degree of gravitational wave from the 
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theory. Later Wilson and Mathews II 211 rediscovered 
one of the approximations of Isenberg and applied 
it to mod el q uasi-equilibria of relativistic binaries (Wil- 
son et al. 13 1, Flanagan [13], Wilson and Mathews |13|). 
Since then the conformal flatness approximation (here- 
after CFA ; also called IWM (Isenberg- Wilson-Mathews) 
theory or CFC (conformal- flatness condition)) has been 
adopted in numerical relativity [H, The simplest 

introduction of CFA is done in (3+l)-decomposition of 
spacetime (see e.g., Gourgoulhon [HI). A general space- 
time metric is written as [15| . 

ds 2 = -a 2 dt 2 + ^ 4 7y (dx l + f3 t dt){dx j + /3 j dt). (1) 

CFA assumes 7y to be a flat spatial metric . Assuming 
the trace of the extrinsic curvature Kij of the spatial 
slice to be zero, we obtain the system of elliptic partial 
differential equations 13|, [l4| : 



pressure. For the comparison of slowly rotating cases we 
use the same equilibria as [17j . For rapidly rotating cases 
we use equilibria computed by COCAL code [l9| assum- 
ing an EOS of polytropic form as p = Kp 1+1 / N , where 
p is rest mass density. It is known that the metric of 
a general stationary and axisymmetric spacetime cannot 
be cast into a conformally flat form on which we perform 
perturbation. Therefore, for internal consistency we use 
equilibria constructed by assuming conformal flatness, 

ds 2 = -a 2 dt 2 + i[j 4 [dr 2 + r 2 d0 2 + r 2 sin 2 6(dtp - Lodtf], 

(5) 

where a, ip and U) are metric coefficients. 

For simplicity we assume that perturbed fluid follows 
the same EOS as in the equilibrium. This suppresses 
g-modes arising from a stratified chemical composition 
and a beta-freezing of nucleons in the perturbed fluid 
El. Thus we have 
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(2) 
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Here Vj means the covariant derivative with respect to 
the flat metric fy. The indices of tensors are raised and 
lowered by f ir k ij is defined as k ij = V l f5 j + V J /3 4 - 
|/ y Vfc/3 . In the source terms we have p H = n a ni,T ab = 
£ (cm^+aW 4 , S* = -b{n c T« - £ a 7rV-^f (ctt*) 2 , 
where n a is the unit normal to the spatial hypersur- 
face of t — const.. Linearized around the equilibrium 
state, these equations are elliptic PDEs for Eulerian per- 
turbations of a, 4> an d ft 3 (j — r,9,ip). For a given 
right hand side, these five equations are formally solved 
for 6a, dip, SP 3 ', by using appropriate Green's functions. 
For scalar variables 5a and Sip, we use Green's function 
1 / 1 f — f 1 1 which is expanded as a sum of scalar spherical 
harmonics. For the vector variable 5/3 J , we introduce a 
vector harmonic expansion of the vector Laplacian oper- 
ator as in Hill [16]. 

The assumptions on the equilibrium stars that we per- 
turb are: (1) stationary and axisymmetric, (2) the equa- 
tion of state (EOS) of the constituent fluid is barotropic. 
For slowly rotating cases, we follow Yoshida and Kojima 
[ijj to construct the equilibria. In their treatment a stel- 
lar structure is computed up to the first order of its rota- 
tional frequency, i.e., physical quantities depend on the 
radial coordinate and they coincide with those of the non- 
rotating ones except for the (t, (^-component of metric 
coefficient [H|. The EOS adopted in 17] is p = Ke 1+1 ' n , 
where e is energy density, K and n are constants, and p is 
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for our slowly rotating cases and 
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for our rapidly rotating cases. A here means Lagrangian 
perturbation of the quantity following it. Relativistic 
enthalpy h is defined as h = (e + p)/p. In the present 
case all the thermodynamic variables e,p, p are expressed 
as functions of h. 

The perturbation equations for fluid variables are the 
equation of rest mass conservation and the equation of 
motion (the conservation of energy- momentum tensor). 
We introduce perturbed fluid variables 5h and <57r a , where 
6 means Eulerian perturbation of a quantity that follows 
it. Momentum density 7r a is defined as ir a = hu a , where 
u a is 4-velocity of the fluid. Then the perturbed rest 
mass conservation is written as 



0. 



(8) 



It should be noticed that V a above is the covariant 
derivative of the four dimensional metric. The perturbed 
equation of motion is expressed as 



5 [^dbiTj - n b djir b ] =0 (j = r, 9, ip) 



(9) 



At a stellar surface where p = 0, Eq.© is replaced 
by a boundary condition. We adopt the conventional 
condition 



Ap = 0, 



(10) 



i.e., Lagrangian perturbation of pressure vanishes. La- 
grangian and Eulerian perturbations are related to each 
other by using the Lie derivative along Lagrangian dis- 
placement vector which is a function of the perturbation 
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of 4- velocity [2l|. Therefore the boundary condition is 
cast into a relation between Sh and Swj. 

The scheme of solutions of eigenmode problem is a 
modified version of Yoshida and Eriguchi [3J . The equa- 
tions for fluid variables (Eq.flSJ), (!§))■ (fTU)) ) are cast into a 
form in which all the fluid variables are on the left hand 
side and the terms containing perturbed metric coeffi- 
cients are on the right hand side. If the right hand sides 
are omitted, the equations in the Cowling approximation 
used in Yoshida and Eriguchi Q are recovered. Once 
the metric perturbation is provided, the right hand sides 
are computed and the system of equations are solved for 
the fluid perturbation and the eigenfrequency. On the 
other hand, Eqs.flU,©,® are inverted when their source 
terms are computed for given fluid perturbations and the 
prescribed gravitational perturbations. Thus the equa- 
tions for fluid and metric variables are iteratively solved. 

RESULTS 



TABLE I. Frequencies of the lowest quadrupole (I = m = 2) 
modes for slowly rotating stars. The EOS is polytropic 
(p — Ke 1+1 ^ n ). a R is the real part of the quasi- normal 
frequency of a spherical star with given compactness M/R 
(taken from Yoshida and Kojima [j?fl )■ o CFA is the eigen- 
frequency computed by the CFA and a 0w is the frequency 
obtained by the Cowling approximation. a' R , o' CFA , ff' are 
rotational corrections to the corresponding frequencies (see 
text for the definition). 



mode M/R *° nFA /* R <x°„>° *' rF Jv' R a' r J< 



f 


0.100 


0.998 


1.26 


0.998 1.02 




0.200 


0.997 


1.15 


1.00 1.01 


pi 


0.100 


1.00 


1.10 


0.994 1.00 




0.200 


0.997 


1.11 


0.991 1.01 




0.100 


1.000 


1.05 


0.997 1.00 




0.200 


0.999 


1.06 


1.00 1.00 



First we compare the eigenfrequencies of low-order 
non-axisymmetric modes of slowly rotating stars ob- 
tained by the CFA, with those by the full GR theory. In 
Table Q] eigenfrequencies computed by the CFA and the 
Cowling approximation are compared with the real part 
of the quasi-normal frequencies in full GR theory. o® CFA 
and cr; are the eigenfrequencies of non-rotating stars for 
the CFA and the Cowling approximation, respectively. 
a R is the real part of the quasi-normal frequency for a 
non-rotating star, a' , , a' , a' are rotational correc- 

° CFA 1 Cw 1 R ■ — -i 

tions to the frequencies defined in [17j ]. 

a = a + ma'Q, (11) 

where f2 is the rotational frequency of the star. It should 
be noted that although the metric of a spherical and 
static background spacetime can be exactly cast into a 
conformally flat form (by introducing the isotropic coor- 
dinate), the perturbed spacetime around it is no longer 
exactly conformally flat. The frequencies computed in 
the CFA are slightly different from the quasi-normal ones 
with all the GR effects taken into account. As expected, 
the errors in frequencies computed in the CFA are smaller 
than the corresponding ones by the Cowling. Although 
the errors in rotational correction (the factor of ft in 
Eg . (TTTT) ) . a', seem slightly larger for the CFA compared 
to the Cowling, the errors in the frequencies for non- 
rotating case (a ) are much larger for Cowling. Thus the 
total error in the frequencies are smaller for the CFA. 
Overall, the error in the frequencies of m = 2 f-, p\-, 
and p2-modes for the CFA is within 1% for the typical 
compactness of a neutron star. 

Next we check errors of cigenmodes due to finite grid 
spacing. Since we work in the finite-difference approx- 
imation of the equations, eigenmodes obtained numer- 
ically depend on the grid spacing. By changing the 



grid spacing, we measure the dependence of the error 
on it. When the finite grid parameter 8 is introduced 
(say, 8 = Ar in the radial grid) , a numerically-computed 
eigenfrequency is developed as a Taylor series around 
the limit of 5 — > 0, thus a [8] — ao + a\5 + ^S 2 + • • •. 
Here 00 is the value of the numerical eigenfrequency com- 
puted by Richardson extrapolation. To evaluate the error 
size we compare the size of each terms in the Taylor se- 
ries. For I = m = 2 f-mode of a non-rotating star with 
M/R = 0.200, we have a = -1.17,01 = 2.18, and a 2 = 
2.41 x 10 -3 , when 8 is the radial grid spacing normalized 
by the radius of the equilibrium star. The rotational cor- 
rection a' defined above is developed in the same way 
and wc obtain a' = 6.71 x lO" 1 ,^ = 4.36 x lO" 1 ,^ = 
-8.20 x 10~ 2 . 

Our code performs (vector) spherical-harmonic expan- 
sion of variables and equations with a truncation of series 
at some value of the degree of the harmonics. For the 
I = m f-mode, a scalar quantity as Sh is expanded by 
{ P m+2k( C0S Q)} ( k = 0, 1, 2, . . . , fc max ) for a given order 
to. Thus the convergence in the angular grid is exam- 
ined through the behavior of the eigenmodes for differ- 
ent A: max . For a rapidly rotating star of polytropic index 
N = 1 and M/ R = 0.1 with its surface polar to equatorial 
axis ratio being 0.65, it is necessary to have fc max > 14 to 
determine an eigenfrequency within 1% of error. A star 
closest to its mass-shedding limit at which the fluid on 
its equator rotates at Keplerian velocity, deforms more 
and the number of harmonic terms necessary to have a 
good approximation of eigenfunctions increases. For a 
star with a smaller rotational frequency a smaller num- 
ber of /c max is sufficient. 

For rapidly rotating cases, we compute sequences of 
the low-order counter-rotating f-modes of the order to = 
2, 3, 4. These modes correspond to I = m f-modes in the 
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FIG. 1. Dimensionless pattern speed uj p (defined as a/m 
where a is the eigenfrequency and m is the azimuthal quan- 
tum number of the eigenmode) of counter-rotating f-modes 
as a function of a stellar rotational parameter T/W. These 
modes become I = m f-modes in the non-rotating limit of 
the star. For a slowly rotating star their pattern rotates 
against the stellar spin. These modes becomes unstable to 
the CFS instability at their zeroes when viscosity is neglected. 
The rightmost points of the mode sequence correspond to the 
mass-shedding limit of the star. The black dots are the neu- 
tral points of the instability computed by SF98. From right 
to left they correspond to m — 2, 3, 4 f-modes. 



non-rotating stars. The pattern speed of these modes are 
retrograde with respect to the stellar rotation. They have 
been of great interest as because they are the most sus- 
ceptible modes to an instability driven by gravitational 
radiation (Chandrasekhar -Friedman-Schutz (CFS) in- 
stability ; Friedman [21], Chandrasekhar [23|, Friedman 
and Schutz [H). 

In Fig[T] the pattern speed u p of counter-rotating f- 
modes with m = 2,3,4 arc plotted for the stellar equi- 
librium sequence of N = 1 with a fixed rest mass. It 
is defined as ui p = a/m. The ratio of the rotational en- 
ergy to the gravitational energy T/W is adopted as a pa- 
rameter characterizing the degree of stellar rotation [25| ■ 
The non-rotating limit of the equilibrium sequence have 
a compactness C (i.e., the ratio of gravitational mass 
to the circumferential radius) value of 0.2. When the 
gravitational radiation is taken into account, an inviscid 
star rotating faster than the zeroes of these modes be- 
comes unstable due to the CFS instability. As is expected 
the difference between a mode sequence computed in the 
CFA and that computed in the Cowling approximation 
becomes smaller for larger m. Introducing gravitational 
perturbations, the pattern speed tends to be larger than 
in the cases with the Cowling approximation. Thus the 
Cowling approximation underestimate the CFS instabil- 
ity It should be also remarked that the slope of the 
pattern-speed curve around T/W = is larger than that 



FIG. 2. Perturbed fluid variables of the counter-rotating 
bar mode near its neutral point of CFS instability (T/W = 
0.0843, the ratio of the polar to the equatorial radius is 0.650). 
The value of functions along each "6 — const." direction 
is plotted as a function of r. Each panel corresponds to 
the following functions ; top-left:<5/t, top-right :<57r r , bottom- 
lcit:Sng/r, bottom-right: STT v /rsind. 




0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4 




.4 0.2 0.4 0.6 0.8 1 1.2 1.4 

R 



0.2 0.4 0.6 0.8 1 1.2 1.4 
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FIG. 3. Metric perturbations for the same model as Fig[2] 
are plotted inside the star. Each panel corresponds to the 
following functions ; top-left:<5a, top-right:^, centre-left:<5/3 r , 
centre-right S/3e/r and bottom:<5/3^/r sin 9. 

around the actual neutral point. Therefore an estimate of 
neutral points by using the slow rotation approximation 
would rather overestimate the instability. 

Eigenfunctions of perturbed fluid and metric variables 
are plotted in Figj2] and [3] respectively (although metric 
variables are computed also in the vacuum domain, they 
are displayed inside the star) . Each curve corresponds to 
the value of a function on a "6 = const." line. Notice 
that the ratio of the polar axis to the equatorial axis of 
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the star is 0.65 and the rightmost point of each curve 
corresponds to an equilibrium surface point. These are 
plots for the counter-rotating bar mode at T/W — 0.0843 
where it is neutrally stable against the bar-mode CFS 
instability. Eigenfunctions are normalized as 5h = 1 at 
the closest surface point to the equator. As a star spins 
up, the eigenfunctions in enthalpy and in the momentum 
change their profiles in such a way that their amplitudes 
are peaked more toward the equator Q. 



CONCLUSION 

We presented a new code that solves the eigenmode 
problem for rapidly rotating stars in GR. The code is 
intended to "go beyond the Cowling approximation" on 
which the precedent eigenmode solvers in GR have been 
mainly developed. Our code adopts the CFA to take 
into account the non-radiative part of relativistic gravity. 
For slowly rotating stars we see that the mode frequency 
computed by the CFA rather improves those obtained by 
Cowling approximation. 

A particularly interesting issue in rapidly rotating rel- 
ativistic stars is their instability. The zeroes of counter- 
rotating modes marks the onset of the CFS instability 
and the slow rotation approximation is rather erroneous 
in evaluating them. Thus the application of the current 
code on this issue is promising. We compared the critical 
T/W available through Fig. 5 in SF98 by reading off the 
numbers there and compared them with our results. For 
the weaker gravity C = 0.1 sequence, the difference of 
the critical T/W for the CFA from that of SF98 is less 
than 1% for to — 3 and 2% for to = 4 mode. For C = 0.2 
the difference is 6% for to — 2 and for to = 4 and less 
than 1% for to = 3. The good agreements for to = 3 
might be coincidental. These results are consistent with 
the analysis done by non-linear simulations 
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